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SUMMARY 


A  study  was  conducted  to  develop  an  implicit  method  for  integrating  the  equations  of 
motion  of  a  lumped-mass  model  of  a  rotor  bearing  system.  The  approach  was,  first,  to  use 
a  Nordsieck-like  numerical  integration  directly  on  the  second-order  equations  of  motion 
and,  second,  to  assume  that  the  forces  and  torques  on  the  rotor  are  functions  of  the  position 
and  velocity  at  the  point  of  application  and  its  nearest  axial  neighbors,  lhis  allows  the 
variables  to  be  arranged  so  that  the  Jacobian  of  the  set  of  nonlinear  equations  is  block 
tridiagonal.  Therefore  the  computational  time  is  proportional  to  the  number  of  elements  in 
the  rotor  dynamics  model  rather  than  to  the  cube  of  the  number.  Numerical  stability  was 
demonstrated  for  any  linearized  homogeneous  mode. 

To  decrease  computational  time,  a  closed-form  solution  to  the  short-bearing  theory 
was  derived  for  a  damper  with  arbitrary  motion.  Explicit  results  were  presented  for  no 
cavitation  and  for  full  cavitation. 

The  vast  amount  of  data  generated  by  the  computer  code  was  displayed  in  a  motion 
picture  showing  an  oblique  view  of  the  rotor  bearing  system.  The  motion  of  the  rotor  could 
be  easily  interpreted. 

An  example  problem  of  a  rotor  accelerating  through  three  critical  speeds  with  19 
mass  elements  in  the  rotor  dynamics  model  took  0.7  second  of  central  processing  unit  time 
per  time  step  on  an  IBM  360-67  computer  in  a  time-sharing  mode.  The  mode  shapes  at  the 
first  and  third  critical  speeds  were  similar  to  the  predicted  mode  shapes  and  occurred  at 
the  predicted  speed.  Because  of  the  unbalance  distribution,  the  second  mode  was  not 
excited.  Above  the  third  critical  speed  the  rotor  bearing  system  operated  as  a  self- 
centering  device.  This  was  also  observed  experimentally. 

The  computer  code,  for  the  first  time,  allows  us  to  look  at  a  complex  rotor  bearing 
system  with  nonlinear  transients  and  displays  the  vast  amount  of  results  in  an  easily 
understood  motion-picture  format.  A  10-minute  16-millimeter,  color,  sound  motion-picture 
supplement  is  available  on  loan. 


INTRODUCTION 

Nonlinear  transients  that  are  important  in  flexible,  rotating  equipment  are  difficult 
to  analyze.  Such  things  as  blade  tip  rubs,  spline  friction,  and  squeeze-film  dampers  are 
difficult  to  predict  with  a  linear  model.  Some  of  the  transients  that  are  important  are 
locked  rotor  starts,  blade  loss,  and  rapid  deceleration  due  to  bearing  failures. 

There  are  two  basic  methods  for  studying  transient  rotor  dynamics.  The  first  me  thou 
is  the  modal  method  (refs.  1  and  2).  It  is  best  suited  to  linear  rotor  bearing  systems  njnning 
at  a  constant  speed.  The  second  method  is  the  direct  integration  of  the  equations  of  mo¬ 
tion.  It  can  be  applied  easily  to  nonlinear  systems  that  are  varying  in  speed.  The  problem 
with  the  direct  method  is  that  it  is  limited  by  either  computer  running  time  or  numerical 
stability. 

The  equations  of  motion  for  rotor  dynamics  cun  be  integrated  directly  in  either  of 
two  ways,  explicit  or  implicit  integration.  The  explicit  integration  method  solves  the 
equations  of  motion  at  the  present  time  for  higher  order  derivatives  and  then  extrapolates 
the  displacements  and  velocities  with  a  Taylor  series  to  the  advanced  time  (ref.  31.  The 
implicit  method  solves  the  equations  of  motion  (implicitly)  at  the  advanced  time  step  for 


the  displacements  and  velocities,  such  that  an  extrapolation  backward  in  time  gives  the 
previous  results. 

The  explicit  method  tends  to  be  unstable  when  the  product  of  the  critical  frequency 
(for  any  mode  numerically  possible)  and  the  time  step  is  large  (ref.  4).  Since  the  highest 
frequency  is  related  to  the  square  of  the  number  of  elements  in  the  rotor  dynamics  moael, 
the  computational  time  will  be  related  to  the  square  of  the  number  of  elements.  Ap¬ 
proximately  five  or  six  elements  seems  to  be  a  practical  limit  to  the  explicit  method  (ref.  2\ 
that  is,  it  can  only  be  applied  to  simple  assemblies. 

In  contrast  to  the  explicit  method,  the  implicit  method  tends  to  be  stable  for  large 
time  steps  (ref.  5);  but  it  requires  the  solution  of  a  large  number  of  nonlinear  simultaneous 
equations  at  each  time  step.  For  every  element  in  the  rotor  dynamics  model  there  are  four 
degrees  of  freedom.  For  each  degree  of  freedom  there  is  an  associated  displacement  and 
velocity.  Therefore  the  total  number  of  nonlinear  equations  to  be  solved  at  each  time  step 
is  eight  times  the  number  of  elements  in  the  rotor  dynamics  model.  The  number  of  com¬ 
putations  necessary  to  solve  these  equations  is  proportional  to  the  cube  of  the  number  of 
equations.  Therefore  the  computing  time  is  proportional  to  the  cube  of  the  number  of 
elements. 

This  study  was  conducted  to  develop  an  implicit  method  for  integrating  the  equations 
of  motion  in  a  reasonable  amount  of  computing  time.  The  approach  is,  first,  to  use  a 
Nordsieck-like  numerical  integration  directly  on  the  second-order  equations  of  motion* 
and,  second,  to  assume  that  the  forces  and  torques  on  the  rotor  are  functions  of  the  position 
and  velocity  of  the  point  where  the  force  or  torque  is  applied  and  its  nearest  axial 
neighbors.  This  allows  the  variables  to  be  arranged  so  that  the  Jacobian  of  the  set  of 
nonlinear  equations  is  block  tridiagonal.  The  computing  time  is  proportional  to  the  number 
of  elements  in  the  rotor  dynamics  model  rather  than  to  the  cube  of  the  number  of  elements. 

besides  the  problems  associated  with  integrating  the  equations  of  motion,  there  is  a 
problem  of  describing  the  nonlinear  damper  force  at  each  instant  of  time  for  an  arbitrary 
orbit.  In  the  past  this  was  done  by  numerically  integrating  the  Reynolds  equation  around  the 
damper  (ref.  6).  This  required  a  considerable  amount  of  computing  time.  As  an  aside,  a 
closed-form  solution  to  the  short -bearing  theory  was  derived  for  a  damper  with  arbitrary 
motion. 


SYMBOLS 

A  coefficients  used  in  partial-fraction  expansion 

a  property  of  shaft  between  mass  stations  defined  in  eq.  (18a) 

b  property  of  shaft  between  mass  stations  defined  in  eq.  (18b) 

C  radial  clearance 

c  property  of  shaft  between  mass  stations  defined  in  eq.  (18c) 

D  diameter 

^This  method  of  numerical  Integration  was  developed  by  Frank  J.  Zeleznik  of  the 
Lewis  Research  Center.  For  a  set  of  fir-'- order  equations,  it  reduces  to  Gear's 
method. 
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modulus  of  elasticity 

force 

torque 

clearance  in  direction  ii 
moment  of  inertia 
index 
index 

axial  length  between  mass  stations 
mass  of  a  rotor  segment 
number  of  rotor  segments 
radial  direction  at  angle  6 
order  of  error  in  Taylor  series 
pressure 

order  of  Taylor  series 
radial  displacement 
stability  matrix 
element  of  S 
time 

time  step 


defined  in  eq.  (5) 


nondimensional  velocity  of  journal  in  rotating  coordinates 


real  part  of  radial  displacement 
imaginary  part  of  radial  displacement 
independent  variable 
axial  coordinate 
given  set  of  constants 
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□  □' 


r 


angle  defined  in  eq.  (43) 
e  eccentricity  ratio 

£  damping  ratio 

9  circumferential  angle 

eigenvalue  of  stability  matrix 
H  viscosity 

u>  frequency 

Subscripts: 


B  bearing 

J  journal 

P  polar 

T  transverse 

+  associated  with  nearest  axial  neighbor  or  root  of  eq.  (45) 

0  start  of  integration 

1  end  of  integration 

Superscripts: 

(')  time  derivative 

( )’  axial  derivative 

(“)  average  or  conjugate 

(k)  kth  time  derivative 

-*■  vector 

unit  vector 


•  i 


NUMERICAL  INTEGRATION 


Given  an  arbitrary  function  Z^U)  whose  derivatives  exist,  t 

expansion  can  be  written: 

h« 


^  r 


L 


lay  lor  i 


zk(t  +  4t)  -  £  lMlizU>(t)  +0q_t 


with  Lagrange's  remainder  of  order  Oq.^.  If  the  arbitrary  function  is  chosen  as 


=  r™ 

k  k. 


the  Taylor  series  for  this  function  becomes 


Zk(t  +  At) 


E  (& 


(t)  +  o. 


where  the  binomial  coefficients  are  definea  as 


(OT 


:'(j  -  "k')T  for  j  4  k 


for  j  <  k 


If  the  form  of  the  remainder  is  chosen  as 


0  -  a,  u 

q  k 


the  Taylor  series  becomes 


Zk(t  +  At) 


2(0 


Zj(t)  +  aRu 


where  is  a  given  set  of  constants  and  u  can  be  aeterminea  from  the  equation  of  motion 
at  the  advanced  time.  The  equation  of  motion  at  the  aavanced  time  is 

£f(i,  r,  if,  t  +  At)  -  0  (7) 


From  the  definition  of  z,  the  various  derivatives  become 


r00  .  z 

(At)k 


Substituting  for  the  various  derivatives  into  the  equation  of  motion  ana  knowing  the 
at  the  previous  time  result  in  the  equation  of  motion  being  a  function  of 

2>(u,  t  +  At)  -  0 


values 

(9) 


This  equation  can  be  solved  for  u  and,  from  this  value  of 
an  error  estimate  to  control  the  time  step. 


u,  the  remainder  can  be  used  as 


N  UMERICAL  STABILITY 


The  analysis  of  the  stability  of  the  numerical  integration  technique  assumes  a  moael 
of  a  rotor  bearing  system  that  is  linearized  at  some  instant  of  time.  The  homogeneous 
equation  of  motion  for  any  mode  is 


••  •  2 

r  +  2oj?r  +  at  r  =  0  (2.0 

where  w  is  the  natural  frequency  and  c  is  the  damping  ratio  for  the  mode.  For  every 
mode  that  is  numerically  possible,  with  nonnegative  damping  ratio,  the  amplitude  must 
either  remain  constant  or  decay  in  time.  The  numerical  integration  is  defined  as  unstable 
if  the  amplitude  grows  in  time. 

From  the  definition  of  z  the  modal  equation  becomes 


2Z2  +  2w  At  ?ZX+  (u  At)2ZQ  =  0  (11) 

Substituting  the  Taylor  series  into  the  modal  equation  at  the  advanced  time  results  in 

q 

.21 


)  -  1)  +  2 1  co  At  t+  (qj  At V 

'  |_2a2  +  20^,0)  At  ?  +  cxQ(u  At)2 

j-o 


(12) 


For  this  value  of  u,  the  Taylor  series  expresses  the  solution  at  the  advanced  time  in  terms 
of  the  solution  at  the  present  time  as 


t'  2cl2  +  2ai“At  a0(<*>  At)2 


Zj  (t)  (13) 


Defining  the  matrix  element  s^j  to  be 


kj 


■(0 


.  j(j  -  1)  +  2jui  At  Z+  (u  At)2] 

2a2  +  2a1wAt  £+  aQ(u  At)2J 


(14) 
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and  the  q-dimension  vector  Z  gives  the  eigenvalue  equation  as 


SZ  -  XZ  (15) 

If  the  |  X  |  >  1 ,  the  amplitude  grows  ano  the  method  is  numerically  unstable.  For  q  =  2, 
the  given  a's  are  ag  -2,  aj  =  3,  and  03  =  1.  In  the  limit  as  iii  it  +  «,  the  maximum 
|x|  -*•  0.  Therefore  if  the  time  step  At  is  much  larger  than  for  a  mode,  the 
amplitude  of  that  mode  will  approach  zero.  If  a  mode  is  to  nave  a  nonzero  amplitude,  wh 
must  be  small.  In  the  limit  as  a  it  +  0,  the  maximum  |x|  1.  Therefore  the  method 

is  numerically  stable  in  the  two  limits. 


EQUATIONS  OF  MOTION 

A  model  of  the  shaft  showing  the  complex  number  representation  of  the  raaial 
displacement  r  is  shown  in  figure  1.  The  radial  displacement  is  the  distance  between  the 
shaft  centerline  and  the  axis  of  rotation.  It  can  be  represented  by 

r  -  x  +  iy  (16a) 

where  the  real  and  imaginary  axes  are  fixed  in  space  perpendicular  to  the  axis  of  rotation. 
The  slope  of  the  shaft  along  the  axis  of  rotation  is 

r’  =  x’  +  iy'  (16b) 

The  position  of  the  shaft  is  then  described  by  r  and  r'  at  all  the  axial  locations. 

The  lumped-mass  moael  of  a  rotor  divides  the  rotor  into  N  segments.  The  mass  and 
inertia  of  each  segment  are  assumed  to  be  concentrated  at  a  point.  These  points  are  then 
assumed  to  be  connected  by  massless  elastic  beams  that  model  the  stiffness  of  the  rotor. 

The  equations  of  motion  for  the  lumped-mass  moael  were  derived  in  reference  7.  The 
sum  of  the  forces  y.F  at  a  point,  must  be  zero,  where 

y.F  =  -mr  +  a_r_  -  (a_  +  a+)r  +  a+r+ 

+  b_r;  +  (b_  -  b+)r '  -  b+r^  +  F  =  0  (17a) 

and  the  sum  of  the  torques  about  a  point  must  be  zero,  where 


y^G  *  -ITr’  +  lulpi’  -  uilpr* 

-  b_r_  +  (b_  -  b+)r  +  b+r+ 

-  c_r^  -  2(c_  +  c+)r'  -  c+rj  +  G  (17b) 
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The  +  or  -  refer  to  the  next  or  previous  axial  location;  and  a,  b,  and  c  are  properties  of 
the  shaft  between  these  locations: 


12EI/L3 

(18a) 

6EI/L2 

(18b) 

>  2EI/L 

(18c) 

If  the  nonlinear  force  F  and  the  nonlinear  torque  G  are  functions  of  displacements 
and  velocities  of  the  point  and  its  nearest  neighbors,  the  and  ^G  are  functions  of  the 
displacements  and  velocities  of  the  point  and  its  nearest  neighoors.  If  the  Taylor  series  of 
the  numerical  integration  technique  is  substituted  into  the  equations  of  motion  for  the 
acceleration,  velocity,  and  displacement,  the  form  of  the  equations  of  motion  becomes 


u,  u+,  u\  u',  up 

(19a) 

5js(u_* 

u»  u+t  ul»  u'»  u+) 

(19b) 

These  equations  form  a  set  of  2N  complex  nonlinear  equations  in  2N  unknowns. 
These  equations  are  solved  by  rewriting  them  as  4N  real  equations  in  4N  real  unknowns 
and  then  using  a  Newton-Raphson  iterative  technique  to  obtain  a  numerical  solution.  The 
Newton-Raphson  technique  assumes  a  solution,  linearizes  the  equations  about  that  solution, 
and  then  solves  the  linear  set  of  equations  for  a  correction  to  the  assumed  solution.  The 
form  of  the  equations  of  motion  results  in  the  linear  set  of  equations  being  block  tri¬ 
diagonal  (fig.  2).  The  block-tridiagonal  form  allows  the  set  of  equations  to  be  solved  in  a 
very  efficient  manner.  The  computing  time  is  proportional  to  N  rather  than  to  N 3  as  in 
the  general  method. 


SQUEEZE-FILM  DAMPER  BEARING 


The  configuration  of  the  squeeze-film  damper  is  shown  in  figure  3.  The  same  con¬ 
figuration  can  be  used  to  analyze  journal  bearings  where  the  journal  and  the  bearir^j  are 
allowed  to  rotate.  If  a>j  is  the  rotational  speed  of  the  journal  and  wg  is  the  rotational 
speed  of  the  bearing,  the  average  rotational  speed  is 


<*>  _  1  (uB  +  tOj)  (20) 

For  a  damper  this  average  rotational  speed  would  be  zero. 

If  C  is  the  radial  clearance  and  r"  is  the  displacement  of  the  journal  center  with 
respect  to  the  bearing  center,  the  clearance  h  in  the  direction  n  is 

h  ■  c  -  r  •  n  (21) 
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If  r  is  the  velocity  of  the  journal, 


If  n  is  at  an  angle  6, 


(22) 


“  k  x  n  (23) 

where  k  is  a  unit  vector  along  the  axis  of  the  damper  bearing  in  the  z-direction.  If  u  is 
defined  as 

J  -  ^Tk  (24) 

then 

w  If  =  (tS  x  r)  •  n  (25) 


The  Reynolds  equation  for  the  short,  plain  damper  journal  bearing  is  presented  in 
reference  6  as 


(h3 

3p\ 

_  -  3h  ^  3h 

dz 

\12U 

%l) 

“  36  +  3t 

If  the  boundary  conditions  in  the  damper  bearing  are 


P(0,  6,  t)  =  0 

p(l,  e,  t)  =  o 


(26) 


(27a) 

(27b) 


and  if  h  is  not  a  function  of  z, 


P 


The  eccentricity  ratio  is 


so  that 


6yz(L  -  z) 

.3 


x 


-  7) 


(28) 


(29) 

(30) 
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If  V  is  defined  as 


e  -  w  x  e 


_  6  LI  Z  ( L  —  Z  )  ~ 

P  =  <■  >  V  .  n  (32) 

CT(1  -  e  •  n) 

A  ~ ► 

The  pressure  is  zero  when  n  is  perpendicular  to  V,  ano  the  pressure  is  greater  than  zero 
when  V  •  n  >  0. 

The  force  on  the  journal  due  to  the  pressure  in  a  segment  of  the  lilm  extenoing  from 
00  to  6j  is 


*--T 


d  /*e1/*L 

e0  *<3 


Pn  dz  d6 


This  expression  for  the  force  can  be  integrated  axially  and  becomes 


3  r&1  (V  •  n)n 

J  (1  —  e  •  n 


(1  -  e  •  n)' 


The  angular  integral  can  be  integrated  by  transforming  the  integral  to  the  complex 
plane  (fig.  4).  Let  V  be  in  the  real  direction,  e  be  at  an  angle  and  n  be  at  an  angle 
6  so  that 

V  =  |v|  (35a) 


e  =  e  e 


(35b) 


(35c) 


Differentiating  the  expression  for  n  yields 

d0  =  -in  ^  dn 

and  using  the  definition  of  the  complex  cosine  yields 


V  •  n 


V(n  +  n  ) 


-*■  *  (ne  +  en 

e  •  n  =  - - - 
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The  expression  for  the  force  becomes 


F  = 


2yDL3Y,  (f  n2(n2  +  l)dn 
k  C 2  /  /  (en2  -  2n  +  e)3 


(39) 


where  the  integral  is  around  the  unit  circle  from  ng  to  nj,  where 

10 

nQ  =  e  (40a) 

i0l 

n1  =  e  (40b) 

The  pressure  is  zero  at  n  =  +  i,  and  the  pressure  is  greater  than  zero  when  /Ce(n)  is  greater 
than  zero. 

For  no  cavitation  the  integral  extends  completely  around  the  journal;  and  by  using  the 
theory  of  residues,  the  force  becomes 


F  = 


(41) 


For  cavitation  the  integral  extends  from  -i  to  +i;  and  by  using  a  partial-fraction  tech¬ 
nique,  the  force  becomes 


F  =  41 


r 

3 

_ 

/uDL3\/v  \ 

r- + Z 

(-n+)^  2  A+i  (~n_)^  2  A_^ 

\ 

l  c2  M 

(1  H-  n2)^1  '  (1  +  n2)^1 

f 

CM 

II 

•n 

_  T  “  _ 

j 

(42) 


where  T  is  defined  as 


T  *  tan 


-i  Vi  -  ui2 

-ReU) 


(43) 


and  T  is  in  the  first  or  second  quadrant.  The  partial-fraction  expansion  coefficients  are 


A 


±3 


n2(n2  +  1) 
(n±  -  nJ*3 


(44a) 
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(44b) 


2n+(2n‘+l)  3A±3 


±2  .  ,3 

(n+  -  n_) 


n±  ~  n? 


6n  +  1 


3A 


±2 


3A 


±3 


±1  .  J  n  -  n_  ,  >2 

(n+  -  n_)  ±  +  (n+  -  n_) 


(44c) 


where  the  roots  of  the  denominator  of  the  force  equation  are 


n+  = 


(45) 


DISCUSSION  OF  EXAMPLE 

The  rotor  bearing  system  described  in  reference  8  was  used  as  the  example  problem. 
This  rotor  bearing  system  consisted  of  a  shaft  with  three  disks  mounted  on  two  axially  pre- 
loaded  ball  bearings  (fig.  5).  The  bearings  were  mounteo  in  a  squeeze-film  damper  journal, 
and  the  journal  had  a  centering  spring. 

The  first  three  critical  speeds  for  the  rotor  bearing  system  without  oil  in  tne  damper 
are  shown  in  figure  6.  All  the  modes  are  bent-shaft  modes.  The  "classical"  hierarchy  only 
applies  to  stiff  shafts;  therefore,  the  classical  mode  shapes  uo  not  characterize  the  actual 
mode  shapes.  The  first  mode,  about  7581  rpm,  classically  would  be  the  cylindrical  mode. 
But  in  this  case,  it  has  a  large  amount  of  bending  outward  near  the  shaft  center.  The  sec¬ 
ond  mode,  about  9235  rpm,  classically  would  be  the  conical  mode.  In  this  case,  it  has  a 
slight  amount  of  bending  outward  near  the  shaft  ends.  The  third  mode,  about  11  248  rpm, 
classically  would  be  the  bending  mode.  In  this  case,  it  has  a  large  amount  of  bending 
throughout  the  shaft. 

Experimentally  the  rotor  was  accelerated  from  rest  through  the  three  critical 
speeds.  The  Lissajous  patterns  for  the  three  disks  were  displayed  on  three  side-by-sioe 
cathode  ray  tubes.  A  motion  picture  was  taken  of  the  CRT's  plus  a  speed  counter.  The 
Lissajous  patterns  at  the  three  critical  speeds  are  shown  on  figures  7  to  9.  The  three 
critical  speeds  occurred  at  about  the  predicted  speeds,  and  the  Lissajous  patterns  corre¬ 
sponded  to  the  three  mode  shapes. 

The  rotor  bearing  system  was  modeled  by  using  19  elements.  The  rotor  was  assumed 
to  have  a  uniform,  in-line  unbalance,  with  a  mass  eccentricity  of  0.00254  centimeter  (1 
mil).  The  equations  of  motion  for  this  system  were  programed  in  FORTRAN  IV  on  an  IBM 
360-67  computer  in  a  time-sharing  mode.  The  equations  of  motion  were  directly  integrated 
by  the  implicit  integration  method,  with  a  fixed  time  step  of  0.12  millisecond.  The  tran¬ 
sient  analyzed  was  the  rotor  accelerated  from  rest  through  the  three  critical  speeds.  The 
rate  of  acceleration  was  8727  rad/sec^.  Each  time  step  took  about  0.7  second  of  CPU  time. 

The  output  at  each  time  step  of  the  calculation  was  displayed  on  a  CRT.  The  display 
showed  an  oblique  view  of  the  rotor  bearing  system,  with  the  bearing  centerline  as  the  ob¬ 
lique  axis.  The  transverse  vibration  is  indicated  by  a  series  of  dots.  Each  dot  represents  a 
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location  of  an  element  in  the  rotor  dynamics  model.  The  scale  of  the  transverse  vibration 
exaggerates  the  amplitude  of  the  vibration.  The  display  on  the  CRT  was  photographed  at 
each  time  step.  These  photographs  were  then  shown  as  a  motion  picture. 

The  computer-generated  displays  on  the  CRT  at  the  first  and  third  critical  speeds  and 
at  a  speed  much  greater  than  the  third  critical  speed  are  shown  in  figures  10  to  12.  The 
mode  shapes  at  the  first  and  third  critical  speeds  were  similar  to  the  predicted  mode  shapes 
and  occurred  at  the  predicted  speed.  Because  of  the  unbalance  distribution,  the  second 
mode  was  not  excited.  The  only  indication  of  the  second  mode  was  a  traveling  wave  super¬ 
imposed  on  the  first  mode  shape.  This  traveling  wave  decayed  when  the  rotor  went  through 
the  third  critical  speed.  Above  the  third  critical  speed,  the  rotor  bearing  system  operated 
as  a  self-centering  device.  The  mass  centerline  coincided  with  the  bearing  centerline. 
Therefore,  the  rotor  displacement  was  uniform  and  in  line,  with  an  amplitude  of  0.00254 
centimeter  ( 1  mil). 

In  conclusion,  this  computer  code  for  the  first  time  allows  us  to  look  at  complex  rotor 
bearing  systems  experiencing  nonlinear  transients  ano  displays  the  vast  amount  of  results  in 
an  easily  understood  motion-picture  format.  A  10-minute,  16-millimeter,  color,  sound 
motion-picture  supplement  is  available,  on  loan,  that  shows  the  test  data  and  the  computer- 
made  motion  picture. 


CONCLUSIONS 

An  implicit  method  for  integrating  the  equations  of  motion  for  a  lumped-mass  model 
of  a  rotor  dynamic  system  was  developed.  The  following  conclusions  were  drawn: 

1.  The  method  was  numerically  stable  for  any  time  step. 

2.  An  error  estimate  was  available  to  control  the  size  of  the  time  step. 

3.  The  computational  time  was  proportional  to  the  number  of  elements  in  the  rotor 
dynamics  model  rather  than  to  the  cube  of  the  number. 

4.  An  example  problem  with  19  mass  elements  in  the  rotor  dynamics  model  took  0.7 
second  of  central  processing  unit  time  per  time  step  an  an  IBM  360-67  computer  in  a 
time-sharing  mode. 

For  the  first  time,  this  code  allows  the  simulation  of  a  complex  rotor  oearing  system 
experiencing  nonlinear  transients  and  displays  the  vast  amount  of  results  in  an  easily 
understood  motion-picture  format. 
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Figure  9.  -  Rotor  passing  through  third  critical  speed. 
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Figure  10.  -  Oblique  view  of  rotor  centerline  as  rotor 
passes  through  first  critical  speed. 
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Figure  11.  -  Oblique  view  of  rotor  centerline  as  rotor 
passes  through  third  critical  speed. 
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Figure  12.  -  Oblique  view  of  rotor  centerline  at  rotor 
speeds  much  greater  than  third  critical. 
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